ON ADAPTIVE EULERIAN LAGRANGIAN METHOD FOR 
LINEAR CONVECTION-DIFFUSION PROBLEMS 



XIAOZHE HU, YOUNGJU LEE, JINCHAO XU, AND CHENSONG ZHANG 

Abstract. In this paper, we consider the adaptive Eulerian-Lagrangian method 
(ELM) for linear convection-diffusion problems. Unlike the classical a posteri- 
ori error estimations, we estimate the temporal error along the characteristics 
and derive a new a posteriori error bound for ELM semi-discretization. With 
the help of this proposed error bound, we are able to show the optimal con- 
vergence rate of ELM for solutions with minimal regularity. Furthermore, 
by combining this error bound with a standard residual-type estimator for the 
spatial error, we obtain a posteriori error estimators for a fully discrete scheme. 
We present numerical tests to demonstrate the efficiency and robustness of our 
adaptive algorithm. 
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1. Introduction 

Convection-diffusion equations oftentimes arise in mathematical models for fluid 
dynamics, environmental modeling, petroleum reservoir simulation, and other ap- 
plications. Among them, the most challenging case for numerical simulation is 
the convection-dominated problems (when diffusion effect is very small compared 
with the convection effect) [43] . Dominated convection phenomena could appear 
in many real world problems; for example, convective heat transport with large 
Peclet numbers [27], simulation of oil extraction from underground reservoirs [23], 
reactive solute transport in porous media |45j . etc. The solutions of these prob- 
lems usually have sharp moving fronts and complex structures; their nearly hy- 
perbolic nature presents serious mathematical and numerical difficulties. Classical 
numerical methods, developed for diffusion-dominated problems, suffer from spu- 
rious oscillations for convection-dominated problems. Many innovative ideas, like 
Upwinding, Method of Characteristics, and Local Discontinuous Galerkin meth- 
ods, have been introduced to handle these numerical difficulties efficiently; see, for 
example, [U [HI HI] and references therein. 

For problems with nearly hyperbolic nature, it is nature to explore the idea of the 
so-called Method of Characteristics PQ; and, this idea has been combined with differ- 
ent spatial discretizations like finite difference (FD), finite element (FE), and finite 
volume (FV) methods. Along this line of research, the Semi-Lagrangian method 
(or, in the finite element context, the Eulerian-Lagrangian method) treats the con- 
vection and capacity terms together to carry out the temporal discretization in 
the Lagrangian coordinate. The Eulerian-Lagrangian method (ELM) gives rise to 
symmetric discrete linear systems, stabilizes the numerical approximation, and the 
corresponding diffusion problems are solved on a fixed mesh [TB, 42J. This method 
and its variants have been successfully applied not only on the linear convection- 
diffusion problem [18; 5U1[5T], but also the incompressible Naiver-Stokes equations 
and viscoelastic flow problems; see, for example, [4l ] 1^ 1 HO I 12] 152 1 154 ] |2"2] . 

Adaptive mesh refinement (AMR) for partial differential equations (PDEs) has 
been the object of intense study for more than three decades. AMR techniques 
have been proved to be successful to deal with multiscale phenomena and to reduce 
the computation work without losing accuracy when solution is not smooth. In 
general, the adaptive algorithm for static problems generates graded meshes and 
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iterations in the following form: 

SOLVE ->■ ESTIMATE ->■ MARK REFINE/COARSEN 

In the ESTIMATE procedure, we usually employ some computable local indicators 
to estimate the local error of the approximate solution we obtain from the SOLVE 
procedure. These indicators only depend on the datum and/or the approximate 
solution, and show in which part(s) of the domain the local error is relatively too big 

or too small in the current mesh. We then MARK these subdomains and REFINE 
or COARSEN them accordingly. 

Local error indicators determine whether the whole adaptive procedure is ef- 
fective or not. Therefore, a reliable and efficient error indicator is the key to a 
successful AMR method and a posteriori error analysis is often used to obtain such 
an error indicator in practice [H [47] . In the context of finite element methods, 
the theory of a posteriori error analysis and adaptive algorithms for linear elliptic 
problem is now rather mature. Convergence and optimality of adaptive meth- 
ods for linear elliptic problems have been proved as the outcome of a sequence of 
work 17, 36, 6j[46j[8]; see the recent review by Nochetto, Sicbert, and Veeser [38] 
and references therein. On the other hand, for the nonlinear and time-dependent 
problems, the theory is still far from satisfactory. A posteriori error analysis for 
nonlinear evolution problems is even more challenging. 

Adaptivity time-stepping is very important for time dependent problems be- 
cause the practical problems sometimes have singularities or are multiscale in time. 
Uniform time step size cannot capture these phenomena. There are considerable 
amount of work in the literature devoted to the development of efficient adaptive 
algorithms for evolution problems. A posteriori error estimators for linear parabolic 
problems are studied in [351 HH HOI EH 131 EH 133] and are also derived for nonlinear 
problems; see [33 [51] 3H1 132] for example. There have been also some efforts for 
extending a posteriori error analysis for the time-dependent Stokes as well as the 
Navier-Stokes equations [51 1311 [317] . In particular, a posteriori error estimates for 
convection-diffusion problems have been discussed in [25[ El 09] . 

It is nature to employ ARM techniques for convection-dominated problems be- 
cause of the complex structures of the solutions and evolution of these structures 
in time. We also notice that spatial mesh adaptivity plays an important role in 
ELM to reduce numerical oscillations and smearing effect when inexact numerical 
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integrations are employed |28j . Early adoption of adaptive characteristic meth- 
ods has been seen since late 1980's [TBI [551 [T5] . A posteriori error estimates for 
characteristic-Galerkin method for convection-dominated problems have been pro- 
posed: Demokowicz and Oden [TB] considered the method of characteristics com- 
bined with Petrov-Galerkin finite clement method for spatial discretization. Hous- 
ton and Siili |26j give an a posteriori error estimator in Z 2 (L 2 )-norm for the linear 
time-dependent convection-diffusion problem using the duality argument. Chen and 
Ji [TT1[T2] give sharp a posteriori error estimations of ELM in L°°(L 1 )-norm for lin- 
ear and nonlinear convection-diffusion problems, respectively. A related L 1 ^ 1 ) a 
posteriori error bound can be found in Chen, Nochetto, and Schmidt |13] for the 
continuous casting problem (convection-dominated and nonlinearly degenerated). 

In the previous error estimators mentioned above, the time residual, which mea- 
sures the difference between the solutions of two successive time steps, is employed 
as a local indicator for temporal error. On the other hand, it is observed that ELM 
is essentially transforming a convection-diffusion problem to a parabolic problem 
along the characteristics. In this paper, we consider a posteriori error estimators 
for ELM, with focus on temporal error estimators. Motivated by Nochetto, Savare, 
and Verdi [37] . in which adaptive time-stepping scheme for an abstract evolution 
equation dtu + $(u) — in Hilbert spaces is analyzed, we can obtain an a poste- 
riori error bound for the temporal error along the characteristics. Combined with 
space adaptivity, adaptive method has been designed and implemented. Our nu- 
merical experiments in Section [5] suggest that the new a posteriori error estimator 
is effective. Moreover, the numerical results also indicate that the proposed error 
estimators are very efficient and can take advantage of the fact that ELM allows 
larger time stepsize. 

The outline of this paper is as follows. In Section [2] we introduce a model prob- 
lem and its Eulerian-Lagrangian discretization. In Section [3j we discuss temporal 
a posteriori and a priori error analysis for the linear convection-diffusion problem. 
In Section [4j we present a posteriori error estimation for the full-discretization. To 
simplify the discussion, we only consider a linear convection-diffusion model prob- 
lem and restrict ourselves to the standard residual-type estimator for the spatial 
error, although the technique discussed in this paper can be potentially extended 
to other problems and spatial error estimators. In Section [5j we consider the im- 
plementation of our adaptive algorithm and present some numerical experiments. 
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Throughout this paper, we will use the following notation. The symbol L 2 
denotes the space of all square integrable functions and its norm is denoted by || • ||. 
Let H k be the standard Sobolev space of the scalar function whose weak derivatives 
up to order k are square integrable, and, let || • \\^ and \ ■ \k denote the standard 
Sobolev norm and its corresponding seminorm on H k , respectively. Furthermore, 
II • ||fe,w and | • \k.ui denote the norm || • \\k and the semi- norm \ ■ \k restricted to the 
domain to C ft, respectively. We also use the notation Hq for the functions that 
belong to H k and their trace vanish on <9ft. The dual space of is denoted by 
H~ k . We use the notation X < (>)Y to denote the existence of a generic constant 
C, which depends only on ft, such that X < (>)CY. 

2. A MODEL PROBLEM AND DISCRETIZATION 

In this paper, we consider the following linear convection- diffusion model problem 
uu 

(2.1) — +b(x,t) -Vu- eAu = f in ft, forie(0,T] 
with the initial condition 

(2.2) u(0) =u in ft 
and the boundary condition 

(2.3) u = on <9ft, forte(0,T]. 

Here ft C R d (d = 1,2,3) is bounded polygonal domain. We assume that b is 
divergence free and vanishes on 9ft for t e (0, T]. We assume that u e L 2 (ft) and 
f eL 2 (Q,T: H- 1 ({I)). 

2.1. The Eulerian Lagrangian method. In this section, we briefly recall the 
construction of the Eulerian-Lagrangian method. As usual, we introduce the char- 
acteristics (particle trajectory) x(t;s, X), where X is the Lagrangian coordinate 
(original labeling) at time s of the particle. It is also referred to the material co- 
ordinate. And x is the Eulerian coordinate at the current time and referred to as 
reference coordinate. Then for the given velocity field b(x,t), the characteristics is 
defind by the following ordinary differential equation: 

dx{t '£' X) =b(x(t;s,X),t), x(s;s,X) = X 
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In the Eulerian coordinates, be chain rule, the material derivative is defined as 
following 

(2.4) Du := du(x(t;s,X),t) = du 

y ' Dt dt dt 



and then we can rewrite (2.1 1 as follow 



(2.5) eAu = f. 

v ' Dt 3 

In order to avoid the deformation of the mesh, ELM uses a fixed spatial mesh 
at each time-step and traces back along the characteristics. The characteristics at 
each time interval [t n -i,t n ] {n — 1,2, . . . , N) has the same original labeling X. 
Assume that x n {t) := x(t;t n ,X) and 

dx n (t) 

(2.6) ^^l = b {x n {t),t), x n (t n )=X, te[t n ^,t n }. 



We discretize the material derivative in (2.5) using the backward Euler method as 
follows 

(2.7) P-D-VM)^^ 

fan 

Here, k n := t n - t n -i, f n := f{t n ), and U n := U n (x n (t n )) = U n {X). It is easy to 
see that x n (t n -i) = x n (t n -i; t n , X) is a function of X in the Lagrangian coordinate. 



We usually refer to (2.7) as the temporal semi-discretization scheme. 



2.2. Feet searching. Our a posteriori error estimate base on the assumption that 
the characteristics are solved exactly, which preserves the determinant of the Ja- 
cobian of the flow. In the numerical analysis, it is difficult to do that and it was 
apparent that a computational realization of preserving the determinant of the Ja- 
cobian for the flow map to be one was crucial. Therefore, we discuss how to integrate 
the following ordinary differential equation for the computation of the character- 
istic feet. The numerical scheme we discuss here has second order accuracy and 
preserves the determinant of the Jacobian of the flow. 



(2.8) -^y(x,t,s) = b(y(x,t,s),s), 



y{x,t,t) 
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where 

(2.9) V-6 = 0. 



The equation (2.8 1 is often called the source-free dynamical systems due to ( 2.9 ). 
For such a system, the solution y{x, t, s) : R d i-» R d is often called the flow map or 
phase flow and it has the following property. 



(2.10) 



det 



dy{x,t,s) 
dx 



= 1, Vie <e 



Vs e 



We shall begin with introducing some popular scheme to compute (2.8) and 



showing that the scheme is indeed volume-preserving scheme for d = 2 but not for 
d = 3. We will then introduce some volume preserving scheme to solve (2.8) for 
d = 3, which is due to Feng and Shang [2~i] . 



In literatures, the following second order numerical scheme for solving (2.8) is 
popular and it seems to first appear in [?]. 



First, we integrate the equation (2.8) using the mid-point rule to obtain : 



(2.11) i (x - y(x, t, a)) = \b (y (x,t, s + Q , a + ^ + 0(k 2 ), 

where k = t — s. 

Second, the right hand side is approximated by a second order accurate extrap- 
olation. Namely, 



(2.12) b(x,s + ^\= \b{x, s) - \b{x, s-k) + 0(k 2 ). 

The following approximation shall also be used : 

- „x / k\ x + y(x,t,s) „,,9, 

(2.13) v\x,t,s + -\= yy 2 ' +0(k 2 ). 

For notational conveniences, let us denote y n = y{x, t, s) and y n+5 = (x + y n )/2. 
Hence we have the following implicit approximations : 



(2.14) 
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To see that the mid-point rule (2.14) results in the volume-preserving scheme, 



let us take the derivative with respect to x for both sides of (2.14). We then obtain 
the following : 



(2.15) 

where 

(2.16) 

and 

(2.17) 



1 



(I-F(s))=A(y n+ *)(I + F(s)) 



F(s) 



dx 



(2.18) 



We solve (2.15) for F(s) and obtain that 

k 



F(s)= [I+-A(y 



Under the assumption that some appropriate finite element space for the velocity 
field is used so that the divergence free condition of the velocity field is imposed in 
the discrete sense, we have 



(2.19) 



tr A(y n+ 2) = 0. 



From (2.19), we conclude that detF(s) = 1 identically. 

The main reason why such an algorithm is popular seems that it has the volume 
preserving property. On the other hand, it is easy to see that the algorithm may 
not result in the volume-preserving scheme for d = 3. Note that for d — 3, under 
the assumption that trA — 0, for H given as follows, 



(2.20) 



H= I 



k 



A 



I- k A 
2 



we have that 

detH = 1 <=> dctA = 0. 
Our purpose here is that by reviewing the volume-preserving scheme in three 
dimension developed by Feng and Shang in |24j , we wish to make sure such a volume 
preserving scheme can be devised in three dimension and confirm our numerical 
scheme can be implemented. As far as the author is concerned, such a special 
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algorithmic detail has not been implemented in the context of the semi-Lagrangian 
scheme. 

The basic idea of constructing the volume preserving scheme for d = 3 is based 
upon the following observation. Following the idea of H. Weyl, we have : 



(2.21) 



b(y(x,t,s),s) 



/ 




dv\ 


\ 


( 


dv 2 

dm 
dv 2 


\ 




OVi 




+ 










dyi 


J 


V 


dy 2 


) 


V 






where 



(2.22) 

(2.23) 
(2.24) 



X3 



db 2 



<'i -J \h(yi,V2,w,s) + —(y 1 ,y 2 ,w,s)) dw 

+ / b 3 (y 1 ,w,x 3 ,s)dw and 
Jy2 



v 2 = 



I-X2 

I b 1 (y 1 ,w,y 3 ,s)dw 



The actual expressions for ^— - , ^— - , and as follows : 

oy 2 dy 3 dy 2 dyi 



(2.25) — - = -b 3 (y 1 ,y 2 ,y 3 ,s) 
dy 2 

dv\ f X2 db\ 

(2.26) — - = b 2 (y 1 ,y 2 ,y 3 ,s) - —±{yi,w,y 3 ,s)dw 
9y3 J V2 dyi 

(2.27) -p- = 6i(yi,ife,ife,«) 



dy 2 

dv 2 
dyi 



(2.2S) p- = - r^( yi ,w, y 3 ,s)dw. 

Jy 2 oyi 
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From this, it is easy to see that (2.21) holds and by construction u 1 and u 2 given 
as follows are divergence free : 



(2.29) 



b 1 = 



(2.30) 



b 2 = 



( ° \ 

dm 

OVi 

V 1 

d H 2 

ov 2 

dyi 

V J 



and 



Let us now denote S\ by the volume preserving scheme for 



ds 



y(x,t,s) 



b l (y(x,t,s),s), 



(2.31) 

(2.32) y(x,t,s) = x 

with i = 1,2, then the following composition is trivially volume preserving : 

(2.33) y n = S% o 

Moreover, assuming Sf is of second order accurate, it is easy to see that the 



composition (2.33) is of second order. The above idea on the composition is from 



Feng and Shang, 



2.3. Full discretization. Discretizing (2.7) with suitable spatial discretizations, 



we obtain fully discrete numerical schemes. In this paper, we focus on the finite 



element method. First, we define the weak forms of (2.5) and (2.7) as usual. We 
define the bilinear form a(-, •) as follow 

(2.34) a{v,w) := e(yv,Vw) \fv,w = Hq(CI). 

We denote the potential by 



1 s I 2 

(w) :— -a(w,w) — - I \Vw\ dx 



(2.35) 

and its Frechet derivative by 5, i.e., 

(2.36) (?(*>), v) := a(w,v). 
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It is easy to see that <p is convex and J satisfies the following inequality 

(2.37) {${w),v-w)<(j>{v)-<t>{w) \/w,v&V. 
In fact, it is well-known that we have the following identity, 

(2.38) a(w, w — v) — (f)(w) — cj){v) + — v\ 2 Vw, «6F, 



where |||d||| := a(v,v). Furthermore, by taking the test function as u — v in (2.40) 



and applying (2.38), we obtain that 

(2.39) (^, u - v) + J>(u) - <f>(y) + \\u - vf = 0. 



The weak forms of (2.5) and (2.7) can be written as 
Du 

(2.40) ( — ,v)+a(u,v) = (f,v) Vt> G V, 

and 

(2 41) { u n -u n -Hx n (t n . 1 )) ^ v) + a([/ „ u) = (/7>) w g ^ 

tin 

respectively. 

On a shape-regular triangulation Th '■— {^i} of 51, we introduce the piecewise 
continuous linear finite element space Vh C V such that 

(2.42) V h := {v h € C(Q) : V h \ n G P 1 ^)^; G Th} H 1/. 

Let 7/j™ and V£ denote the mesh and the finite element space at time t n respectively. 
Then the fully discrete scheme can be written as: Suppose V^~ X G V^ 1 is known, 
find U% G V£ such that 

TT n — TI n ~ 1 (T n (t , ^ 
(2 43) U h V Vn-l)) )Vh)+a(u n i = {f n iV) ^ £ y n_ 

tin 

where G V£ is some suitable approximation of f n = f(t n ). 

3. Error analysis for temporal semi-discretization 
In this section, we focus on the a posteriori error estimation for temporal semi- 



discretization (2.7). For the sake of simplicity, we assume that / = in this section. 
Different from the standard time error indicators for ELM, which measure the 
solution difference along the time direction, our new error indicator measures the 
difference along the characteristics. We will establish a posteriori error estimation 
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based on the new time error indicator and show its efficiency. We will also show 
an optimal priori error bound as a byproduct. Here optimality does not only mean 
the optimal convergence rate which can also be achieved by classic error analysis, 
but also mean the optimal regularity requirement which have not been proved by 
standard a priori error analysis. 



3.L A posteriori error analysis. We first show that the solution of (2.5) satisfies 

Du 



an energy identity along characteristics. In fact, by taking inner product with 



Dt 



on both side of (2.5 1 and applying the following identity 



— 6(u) = (£(«),—), 



we immediately obtain that 



(3.1) 



Du 



Dt 



D 

Dt' 



(u) = 0. 



We can see that the convection diffusion equation preserve the total energy along 
characteristics in continuous level. 

On the other hand, by choosing the test function v — (U n — U n (x n (t n _i)))/k n G 



V in (2.41| and employing (2.38), we have an discrete energy inequality: For any 
integer 1 < n < N, 

(3.2) 



- ^(i-fa-i)) 



0(£/")-^(t7"- 1 (a:"(tn-i))) 



< 0. 



We note that the equality holds in (3.2 ) if there is no temporal discretization error. 
This motivates the following definition: 



(3.3) 



U n -U n - 1 {x n (t rl ^ 1 )) 



k n 



^(t/")-^"- 1 ^"^!))) 



and we can view the computable quantity > as a measure of the deviation of 



numerical solution from satisfying the energy conservation (3.1). We can use £„ as 
a time error indicator for adaptive time stepping in ELM. 

In order to give an a posteriori error bound for our new time error indicator, we 
construct the following linear interpolation 



(3.4) 



t 



U n - 1 (x n (t n ^ 1 )), 
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where 

U n (x n (t n )) = U n (X) and U n - 1 (x n (t n . 1 )) = U n -\x n {t n ^ t n , X)). 

Since the original labeling X of the characteristics does not depend on time, we 
have 



(3.5) 



dU(t) _ U n - [/""^"(in-i)) 



dt 



Vt e (t n -i,t n ] 



Substituting ( |3.5[ ) back to (2.41) and then applying (2.38), we have 

By adding and subtracting U(t), for any v € V, we have 
(3.6) (^, U(t) -v) + <P(U(t)) - <f>(v) + \\\U n - vf = Jft(i), 



where fH(i) is the residual (which does not depend on the choice of test function) 
(3.7) 



m := (d ^fr ' u(t) ~ u ' n) + 0(c/(i)) " 0(cr) - 



We notice that <p is a convexity functional, i.e., 

<KU(t)) < t ^-±<j ) {U n - 1 {x n {t n ^))) + t -^^^U n ) Vt € [*„-!,*„]• 
Hence, $K(i) can be bounded by 

(3.8) m(t) < (t n - 1)^- 

By choosing v = Uit) in ( 2.39[ ), and v = u(t) in (3.6), adding these two equations, 

and using (3.8), we end up with the following inequality 

(3.9) 

~ \\u{x n {t),t) - U{t)\\ 2 + \\u - U(t)f + \\U - C/ n ||| 2 < 2(t„ - t)i n t € (tn-lM 

Then the following upper bound of the new time error indicator holds: 
Theorem 3.1. Let u be the exact solution of (2.1 ) and {Z7"} n=0 be the time semi- 



discrete solution in (2.7). Assume that is defined as (3.3). For any integer 
1 < m < N the folltowing upper bound holds: 



(3.10) 



U(t m ) - U m f + E / " I" - ^"l 2 ^ ^ E ^ 
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Proof. Integration (3.9| in time, we can obtain 



\\u(t n ) - u n \\'+ / \h-u n \l 2 dt < U(z"(<„-i),t„-i) - u n -"{x n {t n ^)) +K£ n 



Since b is divergence free, which implies det Vx n {t) — 1, we have, after changing of 
variables, that 



\\u{t n ) - U n \\ ' + / \\u-U n \\ 2 At< \\u(t n ^)-U n - L \\ +kiU- 



Then (3.101 follows directly from summing up above inequality from n = 1 to m. 
□ 



3.2. An optimal priori error estimation requiring minimal regularity. Tra- 
ditional a priori error analysis for ELM treats the temporal semi-discretization as 
a finite difference method and derive the error estimation based on the Taylor ex- 
pansion (see [IHIIH])- As a result, we obtain an optimal convergence rate but the 
regularity requirement is suboptimal. Taking advantage of the posteriori error es- 
timation of new time error indicator £, we can derive an optimal order priori error 
estimation with minimal regularity requirement on the datum and the solution. 



From the definition of the new a posteriori error estimator (3.3), (2.41) and 



(2.38), we have 



Zn = ^ r t\U n ~U n - 1 (x n (t n _ 1 ) 



By the definition of #(•) and (2.7), we have 



Substitute back into (3.10), we have 



max \\u(t ri 

Kn<N 



U r 



n •' t n-l n= i 
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Set fcmax := max k n and use the elementary inequality 2(v — w, v) > \v\ 2 — \w\ 2 , 

l<n<N 

we obtain 

(3.11) max \\u(t n ) - U n f + V / " \lu - U n fdt 

Kn<N ^ , 

1 N 



Now using the assumption that b is divergence free and summing up (3.11 1 from 
n = 1 to m, we obtain the following optimal a priori error estimation with minimal 
regularity requirement: 



Corollary 3.1. Let u be the solution of (2.1 ) with initial value uq, and {U n } n=1 i 



the time semi-discrete numerical solution of the temporal semi-discretization (2.7 1. 
For any integer 1 < m < N , we have the following priori error estimate 

(3.i2) iko - u m \\ 2 + / III" - u n fdt < -fcLJISMII 2 . 

Remark 3.1 (Regularity Requirement). Only on initial guess. 

4. A POSTERIORI ERROR ESTIMATION FOR FULL DISCRETIZATION 



In this section, a posteriori error estimation for the fully-discrete scheme (2.43) 



is considered. The characteristic feet x n (i n _i) is defined by (2.6| and is computed 
exactly. Besides the new time error indicator, residual type error estimator is used 
as a spatial error indicator. Let R n be the interior residual, i.e., 

(4.1) R n 1= ft - U '^ U h 1 (*"(**-!)) + £A[/ n j on any T e T n 

and J™ be the jump residual also defined, i.e., 

(4.2) := (VEftk - VEftU • 

where e = 8t\ H 8ti is the edge/face shared by t\ €E T% and € 7^ n , v e is the unit 
normal vector of e from n to T2- 

As the analysis of time semi-discretization, we define the linear interpolation 

U h {t) of {Eft}£U as following 

(4-3) U h (t) = t-^K + ^tf r V(*»-i)), 
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Similar with (3.5), we have 



(4.4) 



dt k n 



Now, the following a posteriori error estimation holds 



Theorem 4.1. Letu and {£7£}^L be the solutions of (2.1) and (2.43) respectively. 
For any integer 1 < m < N , there exists a constant C > which does not depend 
on the mesh size and the time step size, such that the following a posteriori error 
bound 



(4.5) 



1 i 171 rt„ 

\\u(t m )-Un\ 2 + ^J2 

n=l"'*n-i 



III lib 171 n £ 

<iKo)-^ o ii 2 + ^fc^ n +c^Mn+E / " \\f(x n (t),t)-m 2 _ a dt 

n=l n=l n=l "'*'»- 1 



holds, where the temporal error indicator and the spatial error indicator r] n are 
defined as 



£n '■— fh ~ 



k n 



(4.6) 



and 



(4.7) Vn :=$>» < := \h\ \\R n f L2 (t) + e £ h e \\ J^ (e 



Proof. For any v £ V and Vh € Vh i we have 



(4.8) - fl,v) + a(Cfl, v) = ( d ^- ft, v-v h ) + a(U^,v v h ). 



Let w = U% in (2.38) and v h = U% in (2.431, we have 



= (fZ- d ^,v-v h )-a(Uj; > v-v h ). 
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Inserting Uh, we have 

- h, u h -v) + 4>{u h ) - 4>{v) + ±fa - uj;\j 2 

= - fh, U h U%) + 4>{U h ) d>(UJt) + (fZ-~,v-v h )- a(UJt, v-v h ). 



Let w = v and v — Uh in (2.381, we have 

+(fh ~ ~, v ~ *h) - a(UZ, v-v h ) + (f^ U h v). 
Then let v = u and v h = lT h l + H%(u - U r h l ) where HJJ : V -> V£ is the so called 



Clement interpolation operator. Moreover, let v = u — Uh in (2.40) and add to 
above inequality, we can obtain that 

~ t \\n-Uh\\ 2 + \lu-Uht + \lu-U h l t 

=(^f K, u h ujt) + cp{u h ) - mi) 

+ (fh-^> u ~ u h - nj*(u - us)) - a(uj;,u - u£ njj( u - ujt)) 

+ (f-fh\u-U h ). 



Then by the definitions ( |4.1[ ), (4.2) and (4.6), we have 
1 d 



2 , 1 in.. „ in 2 , 1 iil, TTnmi 



, u]t h-u h r + -iu-u h \Y + -iu-uz n 

<{t n -t)i n 

+ (R\ [u UJt) - U n h {u U£)) + W - K) - H£(« UZ)]ds 

e Je 

+ (.f-fh,u-Uh) 



Here in last two inequalities, we use the interpolation property of the Clement 
interpolation, the Cauchy Schwartz inequality, and the Young's inequality. The 
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norm || • |_„ is defined by |u||_ a := sup(u, «j)/|ti;fl|. Now, we have 



(4-9) ~ \Hx n {t),t) U h (t)f + - < (tn ~ t)f„ + C Vn + \\\f- fZf-a- 



For any t* G (t m _i,t m ] 5 by integrating (4.9) in time from t m -\ to t rni we have 



-1 i />min(t m ,t* ) 

-\\ U (x(t*),n-u h (t*)\\ 2 + -J t h-m 2 

< l -\\u{x m {t m _ 1 ),t m _ 1 ) - U^- 1 {x m {t m ^))\\ 2 + ~k 2 m U + Ck mVr> 



m.'m(t m ,t* ) 

\\f(x m (t),t)-f?S\\ 2 _ a dt. 



For n < to, by integrating (4.9) in time from i n _i to t n , we have 



I „ , l.o. „. 1 



< 2ii^"(t„-i),«n-i) - ui-\x n {t n _ x ))\v + - 2 <tn + ck nVn + - 1 \\f(x n (t),t) - mi a dt 



tn-1 



Using the assumption that b is divergence free, sum the above two inequalities from 
n = 1 to to, we have 

m r vain(t n ,t*) 
' ' m in 2 



dK»(**),0- W)ll 2 + zE / 11^-^1 

i i m rn t in rt n 

< 2 \H0) - uX + - £ + c£ fc n ^„ + 2 E / !/(*"(*).*) - l-ad*. 

n=l n=l n=l^ t ™-i 



Note that t* is arbitrary, the above inequality implies the upper bound (4.5). □ 



5. Numerical experiments 

In this section, we first introduce the adaptive algorithm we use in numerical 
experiments and then present some numerical results. 

5.1. Adaptive Algorithm. Since the main purpose of this paper is analyze the 
proposed time error estimator, we start with the algorithm for time step size control. 
Let T0L tlme be the total tolerance allowed for the error introduced by the time 
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discretization, according to Theorem |4.1[ that is 

N r t„ 
(5.1) $>^» + 2 (E/ ll/(^ n (i)^)-4"ll-adt) 2 <T0L t 



Therefore, a natural way to satisfy (5.1 ) is to adjust the time step size such that 

tdt • 1 r tn 1 

(5.2) k n ii n <^p and - jf \\f(x n (t), t) - ft\\_ a dt < — V™!^ 

Based on this, now we propose the following algorithm to adjust the time step size 
and control the error under given tolerance T0L time . 

Algorithm 5.1 (H). Given tolerance T0L time , Si € (0,1), 62 > 1, € (0,1), and 
the initial time stepsize Uq. 

(1) Setn:= 1 

(2) Set k n :— fc n _i 

(3) Solve the time semi- discretization problem using time step size k n 



(4) Compute the time error indicator (4.6 1 



(5) Check the error: If (5.2) is not satisfied, k n := 5ik n and goto (2) 

TDL 1 1 

(6) If k n i n < and — jf \\f(x n (t),t) - fZ\\_ a dt < — V0TOL time , 

(7) Let n := n + 1 and goto (2) 

The above algorithm is for adjusting the time step size and controlling the error 
introduce by time discretiziation. When we consider the fully discretization, we 



need to take mesh adaption into account. According to Theorem 4.1 the space 



error indicator (4.7) gives us the elements where the local error is relatively large 
and should be refined. Let T0L space be the tolerance allowed for the error introduced 
by the spatial discretization, then the usual stopping criterion for mesh refinement 



(5.3) k ^Vn < TOL 



space 



As usual, in order to achieve equal distribution of error, (5.3) can be replaced by 
(5.4) Vn < WL f^. 
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The main difficulty for time-dependent problems is the mesh coarsening. We 
need a coarsening error indicator to guide the mesh coarsening procedure. Here, 
we employ the same coarsening error indicator as the one proposed in |10) . Let 
be a coarsening mesh from mesh 7^™, and V£ and V£ be the corresponding finite 
element spaces. It is nature that C Vft. Furthermore, let € Vg and U% be 
the numerical solution at time t n satisfy that 

( U " U ^y i ^-^ , v) + a (U^ v) = (/", v), VveVS 

( h h h K [tn - l)) ,v)+a(UZ,v) = (r,v), V«GVX* 
Subtracting these two equalities and taking v — Ujj — I^U^ £ V^, we have 

Using the elementary equality lab = a 2 + b 2 — (b — a) 2 for the two terms on the left 
hand side of the above equality respectively, we can obtain the following Galerkin 
orthogonal relation 

\m u;:\\ 2 + mi^s - m 2 = \\ur - r n H un 2 + mik - mn 2 

-m-izuzf-Kius-izm 2 - 

Then we have 

Together with the triangular inequality 

IK*n)-^ll<IKO-^ll + ll^-^ll 



and Theorem 4.1 this suggests us the following coarsening error indicator, 
(5.5) Cn := £ G> # := WW - OTI|£ 3(t) + - WiW), 

T 

and the stopping criterion for coarsening 

TOL 



(5.6) Cn < 



coarsen 



T 



This indicator does not depend on the solution Ujj, and allows us to coarsen only 
once after mesh refinement, which avoid the undesirable case that the elements 
which were marked for coarsening must be refined again to reduce the total error. 
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There are several different ways to couple time step size control and mesh adap- 
tion. A simple way is to include space adaptivity in every iteration of Algorithm 
5.1 But this may result in solving a large number of discrete systems and adap- 
tivity may suffer due to the computational overhead. Here we apply a mild way 
to couple the time and space adaptivity which is a modification of the adaptive 
procedure in [33]. At each time step, we first adjust the time step size with the 
old mesh, and then adapt the mesh, and followed by checking the time error again. 
This leads to the following adaptive algorithm for on single time step. 



Algorithm 5.2 (H). Given the tolerance T0L time , TDL space , and T0L coarsen , param- 
eters Si <E (0,1), 62 > 1 and 9 € (0,1). Let and the mesh T^™" 1 have been 
obtained at the previous time step t n -\. 



(1) 



Set Tfo :— 1 j k n :— fc n — 1 and t n — t n _i + k n 



• Solve the fully discrete problem (2.43) on using time step size k„ 

• Compute the error indicators 



(2) While (5.2 1 is not satisfied, 

• Set k n :— 5\k n and t n = t n _i 



k t) 



• Solve the fully discretize problem (2.43) on using time step size k n 

• Compute the error indicators 
End While 



(3) While (5.4) is not satisfied, 



Refine the mesh and produce a refined T£ 



• Solve the fully discretize problem (2.43) on using time step size k„ 

• Compute the error indicators 



While (5.2) is not satisfied, 



Set kyi '. — S\k n and t n — t n ^ H - k^ 



— Solve the fully discretize problem (2.43) on using time step 
size k^i 

— Compute the error indicators 
End While 

End While 



(4) 



Coarsen the mesh to produce a new mesh 7/™ according to (5.6) 



Solve the fully discretize problem (2.43) on T£ using time step size k 
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TDL 1 /**" 1 

(5) Ifk n £ n < and — Jf \\f(x n (t),t) - /f ||_ a d* < — ^TOL tlme , 

• fc„ := <5 2 fc„ 
£rad J/ 

5.2. Numerical results. Now we will present some numerical results to demon- 
strate the efficiency and robustness of our new adaptive method. In our numerical 
experiments, since we use backward Euler to discretize the material derivative, we 
choose 5% = 0.5, 62 = 2 and 9 = 0.5. 

5.2.1. ID Examples. 

Example 5.1. Given the initial condition 

u(x, 0) = e7 



and / = 0, the exact solution for the model problem (2.1 1 (ifb is constant) is given 
by 

U (x,t) = X e -(x-*o-btf/2(A*+2st) 

y ' VA 2 + 2et 

where X is a parameter measures the width of the support of the solution. The 
computational domain is [—1,2], and we give the Dirichlet boundary condition 
u(-l) = u(2) = 0. 



Example 5.2. Given the initial condition 

u(x, 0) = 



1 x < 
x > 0, 



and f = 0, the exact solution for the model problem (2.1 1 (ifb is constant) is given 
by 

uM = \ { erfc i^vS) +exp (?) erfc (I7I) } ' 

2 f°° _ 2 

where erfc(a;) = —= / e s ds is the so-called complementary error function. The 
V 7f Jx 

computational domain is [0, 2], and we give the Dirichlet boundary condition u(0) — 
l,u(2) = 0. 



Figure |5.1| and |5.2| shows the numerical results for Example |5.1| and |5.2| respec- 
tively. We can easily see that the new time error indicator allows larger time step 
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Figure 5.1. Example |5.1[ Left: time step size comparison; Right: 
numerical solution comparison 
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Figure 5.2. Example 5.2 (e = 10 6 ). Left: time step size com- 
parison; Right: numerical solution comparison 
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size, but the numerical solution still maintains the accuracy. This means the the 
new time error indicator is reliable and effective. 

5.2.2. 2D Examples. 

Example 5.3. This example is the benchmark Gaussian-cone problem (cf. |52) ). 
Given the velocity field b = (y, —x) T and the initial condition 



i(x,y, 0) = exp {-[(a; - x Q ) 2 + (y - y ) 2 ]/(2A 2 )} 
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and / = 0. The exact solution is 

A 2 

u(x,y,t) = A2 - ^ cxp{-[£ 2 + f ]/(2A 2 + 4ef)}, 

where x = x — xocos(t) — yosin(t), y — y + xq sin(t) — yocos(f), ^ — owd 

8 

(xo,yo) — ( ,0). The computational domain is [0, 1] x [0, 1] in M 2 . 

The following two tests are generalizations of the one-dimensional problem, Prob- 
lem Ed 



Example 5.4. Given the velocity field b — (1, 0) T , and the initial condition 

u(x 0) - { 1 l f x < - 2 
u\x, y, ) - | Q otherwise, 

and f = 0. T/ien i/ie ea;aci solution is 

u(x, y,t) = \ {erfc + exp (?) erfc } , 

f/ie computational domain is [0, 1] x [0, 1] m K 2 . 

Example 5.5. Given the velocity field b = (1, 1) T 7 and the initial condition as 

, f 1 if x < 0.2 and y < 0.2 

u(x,y,0) = < n . 
v y ; [0 otherwise. 

and f = 0. TTie ezact solution is 

»• " - H rfc (f^D + - © - fc (f£) } &s) + - (f ) - fc (S) } 

and i/ie computational domain is [0, 1] x [0, 1] in R 2 . 



Figure |5.3[ |5.4| and |5.5| show the numerical results for Example |5.3[ |5.4| and |5.5| 
respectively. As shown in the pictures, the new time error indicator allows larger 
time step size and the space error indicator captures the singularity. Overall, the 
adaptive finite element method based on our error estimate is effective and reliable 
for convection dominated diffusion problems. 

6. Conclusion Remarks 

In this paper, we discuss the adaptive ELM for linear convection-diffusion equa- 
tions. We 
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• derive a new temporal error indicator along the characteristics. We are 
able to show optimal convergence rate of temporal semi-discretization with 
minimal regularity of the exact solution. 

• combine the new temporal error indicator with residual-type spatial error 
indicator and obtain a posteriori error estimation for the fully discretization 
ELM. 

• design efficient adaptive algorithm based on the a posteriori error estima- 
tors. Numerical results shows robustness of the new adaptive method and 
allows larger time steps compared with standard temporal error indicator 
for ELM in general. 

For the future work, we are working on generalize the a posteriori error estimation 
for nonlinear convection dominated problems, where the characteristics has to be 
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solved approximately. In this case, ODE solvers that preserves the determinant of 
the Jacobian of the flow will be important. Meanwhile, we will also generalize the 
algorithm to Navier-Stokes equations. 

Acknowledgment. The authors would like to thank Professor Ricardo H. No- 
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